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ABSTRACT 

Fnndamental physical considerations and past tests snggest that there 
may be a problem with discreteness error in N-body methods widely nsed in 
cosmological clnstering stndies. This conld canse problems with accnracy when 
conpled to hydrodynamics codes. We therefore investigate some of the effects 
that discreteness and two-body scattering may have on N-body simnlations 
with “realistic” cosmological initial conditions. 

We nse an identical snbset of particles from the initial conditions for a 
128^ Particle-Mesh (PM) calcnlation as the initial conditions for a variety of 
Particle-Particle-Particle Mesh (P^M) and Tree code rnns. The force softening 
length and particle nnmber in the P^M and Tree code rnns are varied and 
resnlts are compared with those of the PM rnn. In particnlar, we investigate the 
effect of mass resolntion (or eqnivalently the mean interparticle separation) since 
most “high resolntion” codes only have high resolntion in gravitational force, 
not in mass. We show the evolntion of a wide variety of statistical measnres. 

The phase-insensitive two-point statistics, P{k) and i{R) are affected by the 
nnmber of particles when the force resolntion is held constant, and differ in 
different N-body codes with similar parameters and the same initial conditions. 
Phase-sensitive statistics show greater differences. Resnlts converge at the 
mean interparticle separation scale of the lowest mass-resolntion code. As 
more particles are added, bnt the absolnte scale of the force resolntion is held 
constant, the P^M and the Tree rnns agree more and more strongly with each 
other and with the PM rnn which had the same initial conditions, snggesting 
that the time integration is converging. However, they do not particnlarly 
converge to a PM rnn which continned the power law flnctnations to small 
scales. This snggests high particle density is necessary for correct time evolntion, 
since many different resnlts cannot all be correct. Onr resnlts showing the effect 
of the presence or absence of small-scale initial power snggest that leaving it ont 
is a considerable sonrce of error on comoving scales of the missing wavelengths, 
which can be resolved by pntting in a high particle density. 
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Since the codes never agree well on scales below the mean comoving 
interparticle separation, we hnd little justihcation to use results on these scales 
to make quantitative predictions in cosmology. The range of values found for 
some quantities spans 50%, but others, such as the amount of mass in high 
density regions, can be off by a factor of three or more. Our results have strong 
implications for applications such as the density of galaxy halos, early generation 
objects such as QSO absorber clouds, etc. 


Subject headings: cosmology:miscellaneous-gravitation-hydrodynamics- 
methods: numerical-dark matter 


1. Introduction 

A fundamental approximation underlying cosmological N-body simulations is that 
the density field of the universe may be represented by a set of N discrete particles. The 
Universe is thought to be dominated by some form of dark matter whose mass is small 
compared to that of galaxies, probably small compared to that of stars. The formation of 
structure in most scenarios proceeds from a nearly homogeneous mass distribution with 
small perturbations. The formation of stars, galaxies, clusters, and superclusters is probably 
a hierarchical process, built on a spectrum of density perturbations in this background of 
very low mass particles. In this case, the only signihcant discreteness is that which emerges 
as the particles begin to cluster. The particles themselves (atoms, axions, neutrinos, 
WIMPS, or something else) have such a small mass as to be insignihcant compared to that 
of the aggregates being studied in the simulations. Therefore any successful N-body code 
should approximate the solution to the Poisson-Vlasov equations. This can be coupled to 
hydrodynamics, to model the baryonic matter that gave rise to the luminous component. 
In this case, the gravitational force on each simulation particle should be dominated by the 
mean held, with insignihcant discreteness ehects (of which there are many kinds). To trace 
exactly the evolution of the discrete particle distribution may sometimes be in conhict with 
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this approximation. The correct evolution of the density held, sampled by a set of particles, 
is supposed to be obtained by properly adjusting the parameter a which effectively softens 
the gravitational force like oc (r^ + (There are a variety of details of the shape of 

the softening of the potential at short ranges; this is only a representative example). If 
the dark matter is not modeled as collisionless, the dominant gravitational force on the 
hydrodynamic component will also be incorrect. 

The choice of a is a subtle problem unresolved in practice. If the real universe 
consists of galaxies, with no internal degrees of freedom, then a could be simply set to the 
typical size of galaxies ~ 10 kpc provided that the number density of simulation particles 
Usim is comparable to that of galaxies Ugai ~ O.Olh^Mpc”^. For a cosmological volume 
~ (1000/i“^Mpc)^, this requires that the total number of particles is of order 10^, which 
was impossible a decade ago but is quite feasible with current existing computational 
resources. This could appropriately simulate a Universe in which galaxy formation was 
decoupled from large-scale structure, as often assumed in the 1970’s and earlier. Our 
universe does not seem so simple; the prevailing consensus is that our universe is made up 
of dark matter particles, possibly with masses ranging from ttidm = 10“^ to 10^^ eV. If 
this is the case, their number density is 3 x 10^'^(eV/mDM)f^o^^ Mpc“^, where Oq is the 
cosmological density parameter. Then simulation particles do not correspond to physical 
objects but rather provide a method to sample the density held of dark matter in a very 
sparse manner. In this case a must be chosen so as to suppress the artihcial discreteness 
and two-body collision effects of the simulation. These effects are negligible for the real 
universe even if the masses of large stars are typical of the particle size. Therefore they 
should be absolutely suppressed in simulations. Of course the interaction between galaxies 
is important, but galaxy formation is now thought to proceed by a process of hierarchical 
clustering intimately interwoven with large-scale structure formation. When two-body 
effects between galaxies are included, internal degrees of freedom must be included as well. 

This line of thought naturally indicates that a should be close to the mean simulation 
particle separation Igep = since a Igep would manifest undesirable discrete 

sampling due to the sparseness of the density held (Melott 1981, 1990). Hereafter we dehne 
e = an^{^ so that e = 1 corresponds to this choice. While this is consistent with the PM 
simulation methodology with one or more particles per cell (e.g. Hockney & Eastwood 
1988) the choice of e <0.1 which intends to achieve higher spatial resolution is nearly 
universal in P^M and Tree code simulations (e.g., Efstathiou et ah 1985; Suginohara et al. 
1991). So the question is whether or not the choice e = 1 is too conservative and whether or 
not the choice of e <0.1 correctly traces the evolution of the underlying density field of dark 
matter while retaining higher spatial resolution without suffering from discreteness effects. 
This question may have a different answer depending on exactly what analysis is done. We 
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attempt to answer it for the case of some statistical measures. 

The answer to the above question should depend on the specihc problem the simulation 
is addressing. Efstathiou & Eastwood (1981) studied collisionality (one effect of discreteness) 
by means of mass segregation in a P^M code with e <0.1 and found signihcant two-body 
scattering. Peebles et al. (1989) studied this in a PM code and verihed that e ~ 1 was 
necessary to suppress collisionality. Both these tests used reasonable perturbation spectra 
but only had one diagnostic of discreteness effects. It is customary in computational 
physics to use solutions with known characteristics to test codes before applying them to 
new problems. The shock tube is now often being used in this way for hydro-dynamical 
simulations in cosmology, although no one thinks galaxies formed this way. Melott (1990) 
discussed the problem and presented pictorial evidence that lack of mass resolution can 
damage results as much as lack of force resolution. Melott et al. (1997) presented a clear 
case against a use of e < 1 in one dimensional collapse of a plane wave. This study was 
restricted to a formal test case and the precise effect on realistic cosmological simulations is 
not clear. 

Often it has been argued that most fluid elements in hierarchical clustering contract, 
justifying the use of smaller e. In other words, in this point of view, e only need be large 
enough to prevent two body relaxation in collapsed regions. Kuhlman et al. (1996) showed 
that even for regions of relative overdensity 5 > 10 about half the small fluid elements which 
initially contract in one or two directions expand in the third. Thus, due to anisotropy at 
hrst collapse there is plenty of opportunity for scattering due to unphysical discreteness of 
the kind found in Melott et al. (1997). This would decouple results on small scales from 
the initial conditions. 

Another problem has to do with the use of correct initial conditions. Absence of the 
initial perturbations on scales smaller than the particle Nyquist wavelength but larger than 
the force resolution scale results in inaccurate modeling of the merger history of clumps. 
This cannot be compensated by improving the force resolution. The initial power spectrum 
above the particle Nyquist frequency could in principle have a wide variety of forms. 
Therefore, expecting the right answer without putting in the correct initial conditions 
implies that the hnal stage completely forgets that important piece of information; this 
can be only approximately true. As a test of this we generated two PM models identical 
except for the presence or absence of initial perturbations in the range IQkf < k < 64kf, 
where kf represents the fundamental wavenumber corresponding to the simulation box size. 
Comparing the other models to the first one shows the effect of differences in integration of 
initial conditions, while comparison to the latter includes the effect of the presence of this 
part of the initial conditions. 
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The purpose of the present paper is to quantitatively examine the extent to which 
discreteness effects in simulations with e < 1 may change various statistical measures in 
models with a power spectrum not unlike that present in the range of interest of many viable 
cosmological clustering theories. It may be somewhat surprising that such a fundamental 
issue has not yet been examined in detail in the past. The proper comparison between PM, 
P^M and Tree codes (we will generically refer to P^M and Tree codes as High Force Low 
Mass Resolution codes, HFLMR, in this paper) over the dynamical range which we will 
present below has been feasible for almost a decade on supercomputers, and within reach of 
high-end workstations for about hve years. 

An unpublished comparison test coordinated by David Weinberg in the later 1980’s 
attacked some of these issues. It also found general agreement between codes, with 
differences in detail on individual objects and on small scales. There is no contradiction 
between our results and that study. A crucial difference is our exploration of the effect of 
mass resolution while hxing the force resolution. 


2. Models 

We use PM (Hockney & Eastwood 1988; Melott 1981, 1986), AP^M (Couchman 1991), 
and Tree (Suginohara et al. 1991; Suto 1993) codes. The P^M code had adaptive smoothing 
turned off since we wish to compare a standard P^M method. The Tree runs use the hxed 
smoothing length in comoving coordinates, and we set a tolerance parameter 6 = 0.2 which 
is considerably smaller (and thus more accurate) than conventional choices {6 = 0.5 ~ 0.75). 
6 merely controls how far the tree expansion is carried, and thus the accuracy of long-range 
forces. 

The initial power spectrum in all cases was P{k) oc k~^ up to some cutoff, in 
most cases at the Nyquist frequency k = 16/c/ dictated by the runs with the fewest 
particles. Realization of the corresponding density held was generated using the Zel’dovich 
approximation (Zel’dovich 1970) to perturb the particles from their initial lattice 
(Doroshkevich et al. 1980). All the models are evolved in the Einstein - de Sitter universe 
(Do = !)• The comparisons were performed at three different epochs when the nonlinear 
wavenumber k^\ becomes IG/cf, 8A;f, and 4/cf, where /Cni = k^\{A) is dehned as 

r ^nl 

u2(A;„i, A) = A^ P{k)d^k = 1. (1) 

Jo 

In the above A denotes the expansion factor (unity at the initial condition), and these 
values of k^i correspond to the epochs A = 22.36, 42.13, and 92.20, respectively. The latest 
moment we studied corresponds to nonlinearity on the largest scale which does not suffer 
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from finite-volume boundary condition problems (Ryden & Gramman 1991; Kauffmann & 
Melott 1992). The specific runs we tested and the model parameters are shown in Table 1. 
We note that PM codes, which have been extensively used in most physical applications 
with large numbers of particles, are much faster than the other two types. Thus our 128^ 
PM runs took much less CPU time than even the 32^ Tree or P^M runs. The typical 
limitation on PM runs is memory or disk space, while CPU time is the typical limitation 
on P^M or Tree runs. HFLMR codes usually run with e < 1, but we will examine their 
behavior over a range in N and e, pushing them toward e = 1 by increasing the number of 
particles while keeping the absolute scale of force resolution a constant. So far as we know, 
this crucial experiment has not previously been done with HFLMR codes. 

Our primary strategy is therefore to highlight the largely unexplored mass resolution 
issue by varying the number of particles while keeping a constant (it is not exactly constant 
because in fact the shape of the short-range softening function is different in all three codes). 
Within a given code a will be constant, so we can spot trends. Between codes softening will 
be of comparable size. We can dehne rso as the radius where the force drops to 50% of the 
Newtonian value. For the P^M code, rso = 0.92a. For the Tree code = 0.87a. For our 
PM code, rso = 0.95 grid unit, albeit with considerable scatter in the softened zone. 

In most cases we set a = 1 (in units where the box size is 128) in both the Tree and 
P^M codes, and run PM with a 128^ mesh. (The P^M runs had one mesh cell per particle, 
but this is not the factor determining force resolution there.) It is typical of HFLMR codes 
to have a considerably less than Igep , the mean interparticle separation over the entire 
simulation box. In most uses of PM codes, the opposite is true (although in gravitational 
applications it has become customary to have e = 1 or 0.5). 

It is important to note that one cannot represent initial power to higher than the 
Nyquist wavenumber of the particles or the mesh of the FFT used to impose the initial 
conditions, whichever is worse (in fact the latter is rarely a problem in cosmology). 
Therefore most of our runs only have initial power up to kc = IQkf, so that comparisons 
of only the effects of divergent numerical integration can be done. The N = 64^ and 
32^ initial conditions are taken from a subset of the N = 128^ PM runs with the power 
spectrum cut off at kc = IQkf. In order to explore the effect of having initial power at 
higher wave-numbers which is only possible with more particles, we have one PM run with 
kc = bd/cf. We also have two each P^M and Tree runs with a = 0.25 to put some points 
along the pure force resolution axis (which has been emphasized in past tests by most users 
of HFLMR codes). These “extra” runs have 32^ particles, e = 0.0625, implying a = 0.25. 

We now have for the first time a large number of runs using 3 different codes with 
identical initial conditions, identical force resolution, but varying mass resolution. In 
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particular the number of particles varies widely enough to study the same physical system 
with values of e typical of HFLMR codes as well as to study the same system with a PM 
code with similar force resolution and matching mass resolution (e = 1). 

It is important to note time-step limitations. Elementary principles of numerical 
stability require that no particle move more than a fraction of a softening length a, or about 
one-half mesh unit for PM, in a single time-step. This condition was amply enforced for 
all three codes. 


3. Visual Impressions 

The lowest order discreteness effect is merely sampling. In Figure ^ we show a slice 
of the density held of our PM simulation with kc = 16fcf at fcni = 4fcf one cell thick. Both 
pictures show the same conhguration; in one case all of the particles in the slice are available 
and in the other only those from a 32^ subset are used. The discreteness has a major effect 
on the visual impression, corresponding to the noise effect on high-order statistics. A claim 
of hlamentarity in the sparse picture would certainly be greeted with skepticism. Although 
percolation analysis can still detect a signal in such datasets (e.g. Melott et ah 1983) it is 
very noisy. The change in discreteness in redshift surveys is largely responsible for the shift 
in attitude toward superclusters from the 1970’s through the 1980’s. Clumps, hlaments, and 
sheets are progressively harder to see through discreteness noise; many present simulations 
are able to show hlaments but have too few particles to allow sheets to be seen. 

More seriously, if these two different density helds were coupled to hydrodynamics, 
with the dark matter driving the gravitational held, the results would be very diherent. Of 
course, both distributions have the same two-point correlation function, but many other 
things are completely diherent. 

Diherences between the runs described in Table 1 can be seen in Figures ^ to 4 at 
diherent stages of evolution. At a given stage particles which lie within a box at the same 
location are shown. In cases with more particles, only those with the same initial locations 
as those of the 32^ subset are shown and used in calculating all statistics for comparison 
presented hereafter. The following statements are subjective but are made quantitative 
later by cross-correlation studies. 

There is a general resemblance of all plots at the same stage. Since they all have force 
resolution a <1 if this were the only issue we would expect them to be identical down to 
the scale of one tick mark. They are not. For example, the three P^M runs all with a = 1 
(e = 1, 0.5, 0.25 corresponding to N = 128^, 64^, 32^) all look different on scales larger than 


a. Within this series, as the number of particles is increased they come to resemble the 
PM run with = 16/cf, the same initial condition. The same trend is evident in the Tree 
runs. Another observation one can make is that on a scale of 4 cells (the mean interparticle 
separation of the sparsest runs) there are very few apparent differences. The PM run with 
kc = 64fcf looks different from the others, including the PM run with kc = IQkf. The overall 
conhguration is similar, but there are many small-scale features reproduced by none of the 
other models. Consider for example the shape and orientation of the two largest objects 
in Fig the existence of the second object in Fig the partial bridge around (93,36) in 
Fig ^ Power on scales above the force resolution but below the mass resolution of HFLMR 
has made an apparent difference. However it is described, it is evident that the sub-panels 
in each of these Figures are not identical. Again, on scales of the largest interparticle 
separation (4) the differences appear to evaporate. HFLMR codes are based on the idea 
that one can shrink e far below unity to get higher resolution. We see very real differences 
implying integration errors for even the modest e = 0.25. We expect greater differences for 
smaller e. The transfer of power to small scales is strong, but does not make small-scale 
initial power totally irrelevant. 

The reader is referred to the center and bottom center images of Figure 7 in Beacom 
et al. (1991) for another example of the importance of initial power on small scales even 
when it is deep in the nonlinear regime. In that paper, as well as Melott et al. (1993), we 
emphasized the effectiveness of the transfer of power from long to short waves. The general 
position and orientation of objects is determined by initial perturbations on that comoving 
scale and larger, so smaller perturbations are ignorable for this purpose. See also Melott 
et al. (1990), Little et al. (1991), Evrard and Crone (1992). However, here we stress that 
the internal structure of these objects will vary depending on smaller-scale perturbations. 
If we wish to study that internal structure, the smaller-scale initial perturbations must be 
present and properly evolved. 

In all of the papers cited in the previous paragraph, it was concluded that: (a) There 
is a strong transfer of fluctuations from large to smaller comoving scales, (b) There is a 
nearly negligible transfer from small to large scales. This is the reason for the success of 
approaches such as the Truncated Zel’dovich Approximation (Coles et al. 1993) which 
ignore initial small- scale fluctuations, successfully follow large-scale fluctuations into 
the mildly nonlinear regime, producing remarkably accurate results. Although it was 
not emphasized in any of these papers, visual inspection will verify that different small 
comoving-scale perturbations produce different results on these scales. This is easiest to see 
in Beacom et aL(1991) due to the use of a color scale, so that high density regions can be 
inspected (while they are solid black in the typical N-body dot plots). However, the result 
can be verihed by close inspection of any of them. Small-scale power is diluted by mode 
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coupling, but it is not lost. 

In the following sections we quantify some of the differences we found with a variety of 
statistical measures. 


4. Statistical Comparisons 

This section contains the main body of our results. Although we examine a number 
of measures, this must be considered a preliminary investigation. We wish to stress that 
similarity or lack thereof in the studies of one characteristic between various runs should 
not be construed as implying the same thing for some other purpose. The major phase 
differences we hnd later make simple extrapolation based on power amplitudes highly 
questionable. 

Errors are not shown on the quantities presented herein. This is a choice grounded 
in our computational/comparison strategy, which brings out systematic differences. All 
but one of our runs of different code types had identical initial conditions. A 128^ particle 
ensemble (or a 32^ or 64^ subset in those cases where the Nyquist Theorem would not be 
violated) is evolved from this state. Comparisons are done between a 32^ subset which is 
the same in all cases. We have found that using a different subset of 32^ particles, or the 
full ensemble if larger, produces (with one exception) differences so small as to be invisible 
on our plots. Thus our sampling error is essentially zero. It would be possible to create 
some error bars by bootstrap, but they would be meaningless. 

Errors could be created based on cosmic variance-doing mnltiple realizations of each 
box. This, too is essentially meaningless here. If the box size were made larger, the cosmic 
variance would change, but the differences between codes remain (since we are only adding 
linear modes, which almost any code can handle). We could shrink or expand the errors 
at will by this linear procednre, but the systematic differences we hnd are in the nonlinear 
behavior. 

There are two senses in which error is meaningfnl here: 

(1) Do the differences we hnd matter in practice? This depends on the desired accuracy. 
A factor of two was formerly unimportant in many applications in cosmology, but is now 
insnfhcient for many applications. 

(2) What about other initial conditions? The most meaningful kind of comparison 
would be to have runs with other initial power spectra, or non-Gaussian initial conditions. 
This would bring ont how mnch the systematic diherences we hnd vary depending on the 



kind of initial conditions. This must be left for future work. 


4.1. Power Spectra 

We compute our power spectra on a 128^ mesh with CIC cloud assignment of local 
density. In all cases we base our analysis on 32^ of the particles, so that the discreteness 
contribution is exactly the same. We have not subtracted off the discreteness contribution 
from the spectra, because it is not Poissonian at early times and subtraction of such would 
constitute an error. For consistency we also do not subtract it at late times (when it would 
in any case be quite small). 

In Figure ^ we show the spectrum of the initial conditions evaluated on a series of 
progressively hner meshes (with a vertical offset for clarity). The lowest line corresponds 
to the standard choice of showing things up to the particle Nyquist frequency. However, 
it is common practice at very late times to show the autocorrelation function to very 
small radii. This seems to be a contradiction. If the initial spectra had been shown to the 
same resolution, they would look like these (our example is from 32^ particles; with more 
particles the spikes lower in amplitude and move to higher k). Our point is that if the 
code has sufficient force resolution to resolve these spikes (which are not random phase), 
then those are the initial conditions. This study then includes both an examination of 
integration errors as a function of e as well as the effect of the absence of this initial power- 
two independent but practically intertwined effects. 

In Figure ^ are plotted all of the power spectra from the 32^ subset of the runs specified 
in Table 1 at our three output stages. We also show power ratios to help clarify differences. 
Considering the diversity of types of codes and their operating parameters, the agreement 
is encouraging. In the final output stage, major differences are only present beyond about 
half the mesh Nyquist frequency, i.e., k = 32/cf. Some expected trends emerge: the codes 
with the very smallest a have the highest amplitude at high frequencies for a given code 
type. PM codes in general have low amplitude at large k. Other results are unexpected: 
Tree codes also have systematically lower amplitude than P^M codes with the same e at 
large k. For both Tree and P^M codes, for fixed a, as the number of particles increases, 
the nonlinear amplitude increases. This was not anticipated. We speculate that a larger 
number of particles is better able to handle the mode coupling that accompanies nonlinear 
evolution. 

Some support for this idea can be found in the neighborhood of the “kink” around 
k = 16/cf in the earliest stage shown (bottom group) in some evolved spectra in Fig ^ (we 
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temporarily ignore the run with additional high-fc power). All the runs which had the 
same initial conditions fall into two behavior groups here: those with 32^ particles (dotted 
or longdash-dot lines), and those with more particles, without regard to code type, mass 
or force resolution. All runs with 32^ particles overlie one another, and have the lowest 
amplitude in the region of the kink. Recalling that k = 16/cf is the Nyquist frequency of 
the particles, it is entirely reasonable to suppose that this kink is caused by the inability to 
transfer power to higher frequencies limited by the Nyquist Theorem. Note that runs with 
more particles, but no initial power beyond 16fc/ do not show this kink. 

This kink is hidden later; no doubt mode coupling from even longer waves coupled with 
the considerable deformation of the lattice allows this to fill in. But for some codes-those 
with fewer particles-a crucial step in the evolution was handled incorrectly. The fact that 
the power spectrum recovers does not mean the evolution is correct (see also Bouchet et. al 
1985). It probably does mean one can approximate the hnal power spectrum to fairly high 
wave-numbers with poor mass resolution. Note that the point of divergence of the spectra 
in the evolved state is close to the particle Nyquist wavenumber {k = Ibfcf). It is tempting 
to conclude that although spectral agreement is good, one can only have full confidence up 
to this value. Knowing whether this is a general feature or an accident will require further 
study. 


4.2. Two-point Correlation Function 

For a long time the two-point correlation function ^{R) has been a staple of cosmology 
studies in large-scale structure. See Peebles (1980) for an excellent exposition. Formally, 
it is the Fourier transform of the power spectrum discussed in the last section, so no new 
information is present. However, it is packaged differently, and may allow new insights 
through new presentation. 

In Figure we show the evolution of i{R) at our three stages plotted against separation 
R in grid cells (based on 128^ as the box size). ^ is computed (not estimated) by using 
Fourier Transform techniques on a very hne (256^) density mesh. Thus the flattening shown 
at i? = 1 on Figure 7 is due to softening in the dynamical force law, not any smoothing 
due to our method of computing This method is equivalent to using all possible pairs 
(except those at R < 0.5). We show the spikes which exist at early stages of evolution from 
a lattice; they are normally suppressed in publication comparing with data by only showing 
late stages. 

Note that differences at the latest stage are not too great. However, there is a 
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noticeable effect at the second stage near R = 3. All the codes with 32^ particles {P^M or 
Tree) differ from the rest here. We stress that all correlations and spectra are computed 
using CIC weighting (Hockney & Eastwood 1988) from a set of 32^ particles, even when the 
original contained many more. Thus this is a systematic difference and not an artifact of 
sampling. 

In the hnal stage, while there is good agreement at larger radii, for i? < 2 or twice 
the force resolution scale of most of our runs, things begin to break down severely. Both 
dotted lines (e = 0.0625) are high, consistent with the traditional HFLMR emphasis on 
force resolution. However, the P^M code with 128^ particles and e = 1 is also high. We 
hnd evidence that both better mass resolution and better force resolution contribute to 
maintaining the power law to small radii, which provides support for the use of HFLMR 
codes to get these statistics without good mass resolution. However, we do not know 
whether or not it is accidental. Strangely, the P^M run with 64^ particles and e = 0.5 
has the lowest amplitude at small R, occupying a kind of minimum. As we do not have a 
128^ tree run, we cannot be sure of the position of the minimum there. Still, its run with 
e = 0.0625, a = 0.25 overlies the P^M run. At any rate, the non-monotonic behavior of the 
P^M run and the observed effect of both force and mass resolution separately varied as well 
as code differences make it difficult to justify the use of ^ on scales below about 2a. More 
study of these effects is needed. We have no way of knowing whether our results on force 
or mass resolution would continue to change beyond the range in N and a studied here. 
Furthermore, without further study we cannot be sure whether the scale of agreement on 
^ is coupled to the force resolution scale, the mass resolution scale, or some complicated 
combination of the two. We do not consider the differences at large radii to be very 
signihcant. They are of low absolute amplitude, and at the latest stage they are on a large 
enough scale to be suspicious on grounds of boundary conditions (Kauffmann and Melott 
1992). They may reflect some innate differences in the Green function between the codes; 
the Tree code has a very different strategy. 

Conclusions on the two-point correlation: (1) ^{R) appears to agree to within 15% 
down to R ~ 2a, (2) Runs do not converge at smaller separations, (3) It appears that 
HFLMR may approximate effects of high mass resolution reasonably well, insofar as ^{R) 
is concerned, and (4) ^ in N-body simulations for R <2a cannot be regarded as totally 
reliable, since it changes when either A^ or a is varied independently. 
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4.3. Pairwise Velocities 


The statistics which have been discussed up to this point are based purely upon the 
spatial distribution of the mass. To probe any differences in the dynamics we also need to 
consider the velocity held. To compare the evolved particle velocities we compute the mean 
pairwise (PW) velocity and the mean pairwise dispersion (Peebles 1980). We compute the 
mean pairwise velocity using 


< Vu{r) >= 


N{r] 


EE 

* jV* 



( 2 ) 


where the summation is taken over all the pairs with separation r and N{r) the number of 
such pairs. The PW dispersion, crvi 2 i'<")y is similarly dehned by 


(Tv,^{r) 



EE 

* 



( 3 ) 


In Fig H we plot the absolute magnitude of the mean PW velocity and the PW 
dispersion. There are several clear trends. As expected, when e decreases both the PW 
velocity and dispersion increase. This is because the particles can now approach one another 
more closely and scatter more strongly off from one another. This trend is also explained in 
terms of the cosmic virial theorem incorporating the hnite size effect of particles (Suto & 
Jing 1997). There is a more interesting and unexpected trend; as the number of particles 
is increased (keeping the absolute force resolution hxed) the mean PW velocity and PW 
dispersion tends to increase. While the above trend with e might be expected this trend 
which depends upon the mass resolution is surprising and substantiates our claim that 
decreasing the force softening length past the mean particle separation can lead to spurious 
effects. 


4.4. Density Cross—Correlation 


The density cross-correlation 


< 6162 > 

(JiCr2 


( 4 ) 


was introduced by Coles et ah (1993) as a method to help quantify similarity (or lack 
thereof) in the density distribution between various models, where in our case 6 is the 
density contrast on the (32^ or 128^) mesh and a is its rms value. We compute the 
cross-correlation on two different sized meshes. One with 32^ cells which corresponds to the 
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mass resolution scale of our lowest mass resolution runs and another with 128^ cells which 
corresponds to the force resolution scale of our hducial models. This comparison will help 
quantify whether the variety of models tested here produce similar or different small-scale 
structure. 

In Tables 2, 3, and 4 we show the density cross-correlations at the nonlinear stages 
/Cni = 16fcf, 8k{, and 4fcf. The table entries above the diagonal are evaluated on a 128^ mesh 
and the values below the diagonal on a 32^ mesh. In order to assure uniform treatment, all 
density helds are computed from a set of 32^ particles which had the same initial conditions. 
Therefore, any differences are due to systematic differences between the conhgurations in 
the evolved state. We ignore here the additional differences introduced purely by sampling 
effects, as in Figure |^. However, we ran tests to insure that densities computed from full 
ensembles of particles do not give greatly different results, by comparing with full ensembles 
when they exist. 

The cross-correlations on the 32^ mesh are much better than the cross-correlations 
on the 128^ mesh. In fact, most results are close to unity. Thus, the codes agree rather 
well at the mean interparticle separation scale of the coarsest mass resolution runs. This 
coarse-mesh agreement is generally quite high at all stages for all codes with two exceptions: 
(1) The Tree code runs with 32^ particles and (2) the PM runs with a full range of initial 
power kc = 64k{ evolve values of K closer to 0.9 than to 1.0. The hrst may be due to some 
large-scale differences due to the absence of a grid in the Tree code while PM and P^M 
codes share the same grid-based algorithm in computing long-range force. The PM results 
indicate that the absence of correct initial conditions on small scales can still be detected, 
even at late times and on larger scales, when those modes have gone deeply nonlinear. 
These differences, however, are not large. Even the worst agreement on this mesh is better 
than the best agreement on the hne (128^) mesh corresponding to the force resolution of 
most of the runs. 

The hne-mesh cross-correlations (shown above the diagonal in Tables 2-4) display 
a number of interesting patterns: (1) Agreement between codes with the same initial 
conditions, N^ and e monotonically worsens, reflecting amplihcation of differences. (2) This 
worsening is much more rapid between runs with the same small e. For example e = 0.0625 
Tree cross P^M (same e) evolves from 0.95 to 0.63 while for e = 0.25 it goes from 0.96 to 
0.70. For the N = 64^ runs with e = 0.5, it goes from 0.90 to 0.84. This is strong evidence 
that integration errors are being much more strongly amplihed for small e, as suggested by 
Suisalu & Saar (1995), and by Park (1997). As emphasized by these authors, this suggests 
real problems for HFLMR codes. (3) At a given stage, as N and e are increased, runs with 
the same initial conditions converge to one another. PM, Tree, and P^M all agree rather 
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well in this limit. There are also two islands of consistency in Table 4 between the two P^M 
runs with smallest e and a (0.86), and the two Tree runs with smallest e and a (0.87), but 
this agreement does not extent to agreement between two different codes with small e and 
a. The convergence for large N and e does extend to agreement between different codes. (4) 
The run with kc = Q4:k{ is the best standard of comparison for appropriate initial conditions 
in cosmology, since most realistic scenarios have power extending to very small scales, not a 
cutoff. While the comparisons discussed in (1-3) above highlight differences in integration, 
(4) includes the effect of the presence or absence of small-scale power. It is only possible 
to put this power in when a very large N permits it to be sampled by the particles. Time 
evolution shows that cross-correlation of the full-sampled run with other runs becomes 
slowly better, as mode coupling brings down power from larger scales; however, it never 
really achieves a high value with any other run, regardless of the force resolution of that 
run. In fact, the runs with smallest a have the worst agreement with this full-sampled run, 
even though the strategy of small a is to try to follow clustering on small scales. These 
results combined with those on the correlation function/power spectrum appear to indicate 
that this strategy may help produce the correct amplitudes of Fourier components, but 
make the phases worse. Thus it may be that only statistics like the correlation function and 
the power spectrum which contain no phase information show strong agreements between 
the various models. Since the differences on small scales lie in phase information, we will 
look at that in the next section. 


4.5. Phase Correlations 

Many different conhgurations can have the same two-point correlation function (or 
equivalently the same power spectrum) but look completely different. This is due to 
the difference of the phases of the Fourier component of density held which ,^(r) and 
P{k) ignore. Gaussian-random density helds have statistical properties specihed only by 
amplitudes (since random phases have no information), but the properties of non-Gaussian 
distributions are dependent on phases. The end results of N-body simulations and the 
distribution of matter in the Universe are both non-Gaussian, although they may have 
begun from Gaussian perturbations (as our simulation did; see also Suginohara & Suto 
1991; Lahav et al. 1993). 

Testing for phase agreement is another way of checking for agreement between density 
helds - in this case those produced by the N-body code. Given phase information for 
complex coefficients re*® specihed we use 9 = da — Ob io specify the diherence in phase 
between two coefficients. Our measure of phase agreement is < cos 9 >, where the averaging 
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is over spherical shells in k. This measure is 1 for perfect correlation, 0 for uncorrelated 
distributions, and negative for anti-correlated ones. 

Every correlation coefficient in Tables 2-4 now corresponds to a function of k. In view 
of the daunting task of trying to display all these functions, we look instead for patterns. 
Fortunately, strong ones exist. Figure ^ includes all the data for every run cross-correlated 
with every other run at the same stage. All data points are plotted but not joined; many 
overlap. 

They segregate in three groups (two of which begin to approach one another at the last 
stage). (1) The best group (solid squares, left column) corresponds to phase correlations 
among the group with PM, k^ = 16kf, the 128^ and 64^ P^M runs and the 64^ Tree runs with 
the same initial conditions. All runs which £t this dehnition stay in this ” high-agreement” 
group at all stages. This agreement deteriorates with time, albeit slowly. Even at late 
times, < cos6' > is close to 0.5 at the Nyquist Frequency. (2) The worst group (open circles, 
middle column) corresponds to anything compared with the PM, k^ = 64fcf run. Agreement 
with this worsens, then gets better, as evolution proceeds. (3) All other phase comparisons 
lie in a third group (x’s). This includes comparisons between various HFLMR runs, which 
worsen with time, and is shown in the right column. 

This third group (x’s) lies mostly between the other two, but some comparisons are 
outliers, joining the top or bottom group. This is not consistent over evolution but a few 
regularities exist: The outliers (x’s) that join the high correlation group (squares) include 
the two 32^ P^M runs compared and the two 32^ Tree runs compared. Each of these pairs 
are strongly correlated but do not correlate well with anything else, across code types or 
with other runs with the same force resolution, but more particles, within a code type. 
The outlier x’s that join the worst (circles) group are, without exception, all comparisons 
involving one 32^ particle run and one with more particles within a given code, or two 32^ 
runs from different code types. 

These natural groupings correspond to those seen in the density cross-correlations 
in Tables 2-4. Given the close similarity of the power spectrum in many pairs which 
cross-correlate badly, we conclude that while differences in phases and amplitudes both exist, 
the phase differences are primarily responsible for the correlation coefficient differences. 

Statistical measures that are sensitive to phase information will clearly suffer more 
from discreteness that those that are not. The HFLMR strategy may produce a good 
approximation to the power spectrum but not to anything that depends on phases unless 
the density held is smoothed over the interparticle separation. 

To summarize: For given initial conditions, different codes agree well only when 
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the softening approaches the mean interparticle separation. The absence of correct high 
frequency initial conditions (center column) and discreteness errors (right column) prevent 
the general agreement found in the left column. 


4.6. Displacement vs. Density Contrast 


For this series of tests we compare comoving displacement, |Ar|, between the PM 
particles and our standard series of HFLMR runs as a function of the density contrast 
computed using the PM particles. We compute |Ar(5)| comparing to both PM models 
(fee = 16/cf and 64fcf) and at three stages of evolution {k-^x = 16/cf,, 8fcf, and d/cf). 

There are several trends of interest in these plots. First, in agreement with our results 
from the density cross-correlation tests comparison with the PM = 16fcf case appears 
to produce better |Ar| than the PM kc = 64/cf case. Secondly, the trend of increasing 
agreement between the PM and HFLMR codes as e is increased also appears to hold here 
as well. Thus, we see that increasing the number of particles present in the HFLMR codes 
without changing the absolute force resolution causes their time integration to converge 
more and more toward the PM code result. In Figure 12, the best overall agreement with 
the PM code is the P^M code based on 128^ particles. 


The differences that do exist between different codes with the same initial conditions 
are amplihed as the the particles continue to move apart from one another. Thus, comparing 
the top three panels of Fig Fig and Fig|^, we see that amplihcation of errors is more 
rapid for smaller e, even when the force resolution a is held constant. In Figures |n, 
and 1^ the comparison is done against displacements of particles compared with the PM kc 
= 64fcf run. In this case they are uniformly bad. Although the integration errors are greatly 
reduced by having more particles, the absence of the extra power in the initial conditions 
prevents convergence. None of the runs, regardless of its parameters, comes close. 


4.7. Cluster Analysis and Genus 

In cluster analysis, percolation properties are used to examine the properties of 
connected regions. Criteria of density thresholds or the close approach of particles are 
used to decide which clumps are connected. This has been used to study the formation of 
structure in simulations and in redshift surveys. Percolation properties strongly depend on 
sampling: the more particles in the model, the easier percolation. Obviously percolation in 
the model corresponding to Fig. lb is easier than in Fig. la. 
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There are less obvious differences in the number and other properties of overdense 
regions. Fig. shows three statistics calculated for the subsamples with equal number of 
particles {N = 32^). The left panels show the fraction of mass in the regions having densities 
above a given density threshold for all models at three stages of evolution. The middle 
panels show the fraction of volume occupied by these regions and the right panels show the 
total number of overdense regions as functions of the density threshold. The clumps were 
identified by using the percolation code described in Yess & Shandarin (1996). We show the 
high density part of these distributions. The mass fraction runs from about 0.1% to 10% 
and the number of clusters from 3 to 300 - 1000. There are significant differences between 
the models. Unfortunately, the patterns are different at every stage therefore it is difficult 
to describe them. The differences can be seen more clearly when the ratios are plotted. To 
do so we choose PM (fee = 64fcf) as a fiducial model and plot the logarithms of the ratios in 
Fig. ig. We also extend the range of densities toward lower densities in this figure so the 
convergence of the models can be seen. 

At the stage k^i = 16/cf all models converge at about p < 30. In the range 30 < p < 100 
the fiducial PM64 model has more clumps and more mass in clumps than any other model, 
however at the very high densities p > 100 it has fewer clumps and less mass in the clumps 
than any other model. The typical differences are 2 -3 times. One could expect similar 
errors on small scales when studying the early stages of clustering, such as Lyman-a 
absorbers or high-redshift galaxy formation, or the inner parts of dark matter halos. 

At the stage ka\ = 8fcf the differences between models are the least of all three 
stages. All models except the Tree code with N = 32^, a = l,e = 0.25 have more mass 
in clumps; they are also more numerous and occupy greater volume. The P^M model 
with N = 128^, a = 1 has the highest values of all statistics and the Tree code with 
N = 32^, a = 1, e = 0.25 has the lowest. The P^M and Tree codes with e = 0.0625, a = 0.25 
are significantly lower than P^M model with N = 128^, a = 1. 

This pattern generally persists at the stage fcni = 4/cf, however the amplitude of the 
differences is greater. The P^M model with N = 128^, a = 1 is again among the highest but 
the P^M and Tree codes with best force resolution e = 0.0625, a = 0.25 are now comparable 
to it. All the models converge at about p < 100. When the densities are computed on the 
32^ mesh the differences are considerably smaller but still are noticeable at the fcni = 16fcf 
stage in Fig. However, they almost vanish at the final stage. This variation of peak 
density with particle number is completely consistent with results presented by Craig (1997) 
on the central density of halos in N-body simulations. After this paper was submitted, 
a preprint appeared (Moore et al. 1997) which shows the same effect in CDM models at 
higher density levels. 
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The genus is a measure of the connectivity of structures widely used in cosmology. It is 
specihcally highly sensitive to phase information. All Gaussian distributions have the same 
distribution of genus as a function of volume fraction, with the amplitude merely being 
determined by an integral over the power spectrum. We show the genus here computed 
using the method of Weinberg et. al (1987). The genus plotted for the last stage of the 
evolution in Fig. ^ is essentially in agreement with the other results: there are differences 
up to 15% between different codes on the scale of the force resolution (right panel) and 
almost no differences on the scale of the mean particle separation (left panel). In fact with 
this smoothing the N-body results are shown to reproduce well the quasi-linear prediction 
for the genus curve. See Matsubara (1994) and Matsubara & Suto (1995) for details. 


Cluster analysis has shown that if the clumps were selected using isodensity surfaces 
calculated on the force resolution scale then different codes would predict very different 
numbers of clumps, their total volume and the total mass in the clumps. Differences as 
large as factor of 2 or 3 are typical for the few hundred largest clusters. The vertical scales 
in Fig. 1^ are about 2 -3 times greater than the horizontal scales depending on the statistics 
and stage. This means that if one reverses the question and asks what is the difference 
in the density thresholds provided that same number of clumps (or total mass) selected 
then the difference in densities would be about 30 - 50%. If the densities are computed 
on the scale of the mean particle separation then all codes agree with each other at stage 
fcni = 4fcf with about 20% accuracy in terms of the numbers, volumes and masses at given 
density thresholds or with accuracy of 5% in terms of density. These differences are based 
on clustering, not sampling, since we used the same number of particles for each. However, 
one can get an idea of the size of a possible sampling effect by comparing the results shown 
here for PM with a fully sampled PM density held using all 128^ particles. The differences 
in these statistics between sparse and full PM samples are much smaller than those we 
hnd between different codes. Due to the non-monotonic nature of the changes we see with 
evolution, we cannot make any predictions about how various codes and strategies would 
compare at later stages of evolution. 


5. Discussion 

A number of the statistical results presented here show only moderate differences 
between the various codes and our choices of a and Igep- This is particularly true for the 
two-point statistics, P{k) and ^{r), which contain no phase information. On the other 
hand, for | < Vi 2 (r) > | and (^^^^(r) the dependence on a and Igep is weak, while the 
dependence is strong for the phase correlations, the mass density and number of clumps. 
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as well as the phase correlations and the density cross-correlations. This implies that the 
distribution of the phases may be the pivotal piece of information that is needed to resolve 
this issue. Unfortunately, no satisfactory measure of the distribution of the phases exists. 

We cannot easily compare our results with others, because no study has included 
the variety of codes, mass resolution, and inclusion of normally absent small-scale 
perturbations. The most similar work is the unpublished study coordinated by D.H. 
Weinberg. An analysis similar to our Figure 8 produced very similar results. Efstathiou et 
ah (1988) studied evolution of power-law initial conditions in a P^M code, using primarily 
low-order or averaged statistics. They did hnd fluctuations of order 50% in the rescaled 
multiplicity function (their Figure 9); while not equivalent, this is roughly compatible with 
our results in the central column of Figure 16. It is interesting to keep in mind that in 
plasma physics (e.g., Hockney & Eastwood 1988) and in beam dynamics (Habib 1997) it is 
widely recognized that Isep a is a requirement to accurately model a collisionless system. 
This is in strong contradiction to usual practice in cosmological simulations where Isep 3> a. 
This is a survival in methodology from the original cosmological N-body simulations which 
assumed to evolve the distribution of galaxies and not the dark matter distribution (in 
fact, they were generally called ’’galaxy clustering simulations” at the time). Note that 
we cannot appeal to which code produces the best agreement with observations to settle 
our dilemma. On small scales, where we hnd signihcant differences, there is an unknown 
correspondence between mass and light in the Universe. This issue must be settled by 
hrst principles, i.e., agreement on which equations are to be solved and the most accurate 
methods to solve them. 

There is one crucial piece of evidence in our study. No one would dispute that the 
performance of any of our codes would become more accurate as more particles are added, 
if the force resolution is held constant. As is shown in many of our tests, this causes both 
Tree and P^M results (given the same initial conditions) to converge to our PM run. For 
low-order statistics, such as the power spectrum and correlation function, a small softening 
parameter appears to replicate the effect of more particles (Although the results are not 
identical so it not possible to determine how small the softening parameter can be made 
without introducing artihcial collisionality into the simulation). But on the other hand, it 
is clear that decreasing the softening parameter has profound effects on any measurement 
which depends upon the phase information. 

The sources of the observed errors are not random but systematic in origin. It is 
important to stress that there are two numerical issues; (a) Two codes with identical initial 
conditions and similar a and Isep should produce similar statistical results for all statistical 
measures of the evolved particle distribution and not just a subset. Our results show a 
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convergence in the time integration between various codes as a ~ Igep- Unfortunately it 
is not standard in cosmology (except in about half the PM applications) to use a ~ Igep- 
(b) The trouble is that the “correct” initial conditions correspond to our case kc = 64kf. 
Thus, it appears that the HFLMR codes are not able to correctly model the kc = Q4:k{ 
solution despite the past claims in the literature and will produce solutions which only 
agree to Igep and not down to the force softening length a. This, of course, creates a real 
dilemma when it comes to simulating truly hierarchical structure. PM codes are fast, but 
the number of particles are limited by disk storage requirements. So-called adaptive mesh 
rehnement algorithms will suffer the same fate as the HFLMR methods tested here since 
they are unable to input the true power spectrum down to arbitrary scales initially. This 
leaves us with nested-grid PM methods (Villumsen 1988; Anninos, et al 1994; Splinter 
1996) or related strategies (Moore et al 1997) as a surviving option. These methods allow 
for both high force and mass resolution simultaneously, thus guaranteeing that a ~ Isep is 
satished on the smallest scales and that the true initial power spectrum is input ab initio. 
Local adaptive mesh rehnement schemes may offer some advantages if the rehnement is 
carefully designed and tested. This may allow some reduction of the integration error we 
have demonstrated, but cannot be a substitute for inclusion of the proper initial conditions 
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Fig. 1.— We show density fields constructed from a slice of an evolved state of our PM 
simulation. The parameters of the two images are identical, except that (a) is constructed 
using only those in the 32^ subset (see text) of particles which happen to lie within the 
slice. Image (b) is from all of the 128^ particles which he within the slice. Although the 
inferred two-point correlation would be identical, the visual impressions are very different, 
as would be the noise in higher-order measures of structure. The general awareness of 
superclusters which arose in the 80’s is primarily the result of an analogous change in the 
signal-to-discreteness noise ratio in redshift surveys. 
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Fig. 5.— The power spectrum constructed using the 32^ particles of the initial conditions 
for all except our = 64A;f runs, as evaluated on mesh of size 32^ (solid), 64^ (dots), 128^ 
(shortdash), and 256^ (longdash), shown with vertical offset. The normalization is such that 
a Poisson distribution of points of the same number of particles as mesh cells on each mesh 
would converge to P{k) = 1. Normally, such spectra are only shown up to the particle 
Nyquist frequency, as in the solid line. The spikes are a result of the lattice of particles 
which is deformed to provide the initial conditions, and are of course not random phase. In 
HFLMR codes, this part of the initial conditions is evolved, with nnknown conseqnences. 
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Fig. 6.— Right panel: The power spectrum of three successive evolved stages for all 
our simulations, evaluated from 32^ particles on a 128^ mesh (their force resolution). The 
normalization is such that a Poisson distribution of 128^ particles would converge to P = 1 
Left panel: The ratio of the power in a given model to that in our hducial kc = 64/cf PM run 
at the hrst (bottom) and last (top) evolved stage. The light solid line is the kc = 16/cf PM 
run; the heavy solid line is the kc = Q4:kf PM run. Other heavy lines are P^M runs and light 
lines are tree code runs. Dotted lines: e = 0.0625 runs with 32^ particles. Longdash-dot: 
e = 0.25 runs with 32^ particles. Shortdash: e = 0.5 runs with 64^ particles. Longdash: 
e = 1.0 run (P^M only) with 128^ particles. 
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Fig. 7.— The two-point correlation function of all the mass in the simulations is shown 
for the same three stages shown in the power spectrum plot. Line types here and in all 
subsequent plots match Fig. 6. The spikes which exist at early stages are a relic of the 
lattice upon which the particles began, which is insufficiently deformed to be suppressed 
then. The last stage corresponds to nonlinearity sufficiently strong that boundary conditions 
would become a problem if the simulation were continued further. 
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Fig. 8.— The pairwise velocity and pairwise dispersion, as defined in the text in Section 
4.4 are plotted for each stage, using the usual line types. The force resolution scale is 1 for 
most of our runs, 0.25 for two. The mean interparticle separation of the sparsest runs is 4; 
convergence is not reached until approximately this scale (rather than the force resolution 
scale). 
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Fig. 9.— The successive plots contain all the data for the averaged phase agreement between 
all of the simulations runs at the same stage; < cos6 > is dehned in the text and is 1 
for agreement and 0 for uncorrelated phases. Individual functions are not distinguishable 
here, but it is clear that they fall into classes shown here. The left column corresponds to 
all runs with good mass and force resolution and the same initial conditions. The center 
column corresponds to anything cross-correlated with the run which continued the power-law 
perturbations to wave-numbers impossible for the HFLMR codes. The right column contains 
everything else. The high values found in the right column correspond to two HFLMR P^M 
runs compared and two HFLMR Tree runs compared-which agree within the pair but not 
with anything else. 












Fig. 10.— The scatter plots contain the difference in positions between particles which 
had the same initial position and velocity in the PM run with kc = 16/cf and positions in 
other runs with identical initial conditions (and in the PM, kc = 64/cf run, upper right); 
The distance is plotted against the local density of the simulation in the fiducial PM run. 
Shown here for the state /cni = lO/cf. The line and error bars represnt the mean and standard 
deviation of the distance; the mean is sometimes not centered due to superposition of dots. 
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Fig. 16.— Three rows show three stages in the evolution of the models. The left column 
shows the fraction of mass, the middle column - the fraction of volume, and the right column 
- the number of disjoint regions having densities greater than the density threshold shown 
on the horizontal. The densities were calculated on the 128^ mesh. 
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Fig. 17.— The ratios of the same parameters plotted in Fig. 16 to the hducial PM model 
with N = 128^, and kc = 64. 
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Fig. 19.— The genus is plotted for the hnal stage kni = 4: The left panel shows genus of the 
density helds smoothed on the scale of the mean particle separation; the right panel shows 
genus for the helds smoothed on half of the force resolution scale (only the part corresponding 
to high densities is shown). 
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Table 1. Model Parameters for the Test Cases 


Code 


N 

l-sep^ 

^force 

PM 

kc = 16 

128^ 

1.0 

1.0 


kc = 64 

128^ 

1.0 

1.0 

P^M 

kc = 16 

32^ 

4.0 

0.0625 


kc = 16 

32^ 

4.0 

0.25 


fcc - 16 

64^ 

2.0 

0.5 


kc = 16 

128^ 

1.0 

1.0 

Tree 

kc = 16 

32^ 

4.0 

0.0625 


kc = 16 

32^ 

4.0 

0.25 


kc = 16 

64^ 

2.0 

0.5 


^Isep is the mean particle separation in grid cell units. 

^^force i® i^ units of the mean particle separation, Isep- Note that for most of the runs a = ^ ^sep = only for two does a = 0.25. 




Table 2. 

Cross- 

-Correlations at kni 

= 16 






PM 


P^M 




Tree 


kc = 16 

kc = 64 

128^6 = 1.0 

64^e = 0.5 e 

- 0.25 € 

- 0.0625 

64^6 = 0.5 

€ - 0.25 

e - 0.0625 

PM (kc = 16) 

_ 

0.31 

0.94 

0.93 

0.77 

0.76 

0.97 

0.74 

0.74 

PM (kc = 64) 

0.92 

— 

0.31 

0.30 

0.26 

0.26 

0.31 

0.26 

0.26 

P^M 128^(e = 1.0) 

1.00 

0.91 

— 

0.95 

0.72 

0.72 

0.95 

0.69 

0.69 

P^M 64^(6 = 0.5) 

1.00 

0.91 

1.00 

— 

0.81 

0.81 

0.97 

0.78 

0.77 

P^M(e = 0.25) 

0.99 

0.91 

0.98 

0.99 

— 

0.98 

0.81 

0.96 

0.96 

P^M^ = 0.0625) 

0.98 

0.91 

0.98 

0.99 

1.00 

— 

0.81 

0.95 

0.95 

Tree 64^ (e = 0.5) 

1.00 

0.92 

1.00 

1.00 

0.99 

0.99 

— 

0.79 

0.79 

Tree(e = 0.25) 

0.98 

0.90 

0.98 

0.98 

1.00 

1.00 

0.99 

— 

0.99 

Tree(e = 0.0625) 

0.98 

0.90 

0.98 

0.98 

1.00 

1.00 

0.99 

1.00 

— 




Table 3 

. Cross- 

-Correlations at kni 

= 8 






PM 


P^M 




Tree 


kc = 16 

kc = 64 

128^6 = 1.0 

64^e = 0.5 e 

- 0.25 € - 

- 0.0625 

64^e = 0.5 

€ - 0.25 

e - 0.0625 

PM (kc = 16) 

_ 

0.48 

0.87 

0.84 

0.60 

0.58 

0.91 

0.55 

0.56 

PM (kc = 64) 

0.94 

— 

0.47 

0.46 

0.40 

0.40 

0.46 

0.37 

0.38 

P^M 128^(e = 1.0) 

1.00 

0.94 

— 

0.88 

0.59 

0.57 

0.86 

0.53 

0.54 

P^M = 0.5) 

1.00 

0.94 

1.00 

— 

0.65 

0.63 

0.90 

0.58 

0.59 

P^M(e = 0.25) 

0.97 

0.92 

0.97 

0.98 

— 

0.92 

0.64 

0.86 

0.87 

P^M^ = 0.0625) 

0.96 

0.92 

0.97 

0.97 

1.00 

— 

0.62 

0.83 

0.84 

Tree 64^ (e = 0.5) 

1.00 

0.94 

0.99 

1.00 

0.98 

0.98 

— 

0.61 

0.61 

Tree(e = 0.25) 

0.96 

0.92 

0.96 

0.96 

0.99 

0.99 

0.97 

— 

0.93 

Tree(e = 0.0625) 

0.96 

0.92 

0.96 

0.96 

0.99 

0.99 

0.97 

1.00 

— 
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Table 4. Cross-Correlations at kni = 4 



PM 


P^M 




Tree 


kc = 16 

kc = 64 

128^6 = 1.0 

64^e = 0.5 e 

- 0.25 

€ - 0.0625 

64^6 = 0.5 

€ - 0.25 

e - 0.0625 

PM {kc = 16) 

_ 

0.65 

0.83 

0.80 

0.64 

0.61 

0.84 

0.52 

0.49 

PM {kc = 64) 

0.96 

— 

0.63 

0.63 

0.57 

0.55 

0.63 

0.47 

0.44 

P^M 128^(e = 1.0) 

0.99 

0.95 

— 

0.89 

0.66 

0.64 

0.84 

0.50 

0.48 

P^M 64^(6 = 0.5) 

0.99 

0.96 

0.99 

— 

0.68 

0.66 

0.84 

0.49 

0.47 

P^M(e = 0.25) 

0.96 

0.94 

0.97 

0.97 

— 

0.86 

0.67 

0.70 

0.67 

P^M(e = 0.0625) 

0.96 

0.94 

1.00 

1.00 

1.00 

— 

0.63 

0.63 

0.63 

Tree 64^ (e = 0.5) 

1.00 

0.96 

0.99 

0.99 

0.97 

0.97 

— 

0.60 

0.56 

Tree(e = 0.25) 

0.92 

0.90 

0.92 

0.92 

0.96 

0.96 

0.94 

— 

0.87 

Tree(e = 0.0625) 

0.92 

0.90 

0.91 

0.91 

0.96 

0.96 

0.94 

1.00 

— 














